## Requires two arguments: ctrl.igv file and exp.igv
#outputs normalized igv file

import sys, numpy
from math import *
from numpy import *

TActrl = 0; Readsctrl = 0; HitSitesctrl = 0
inputreads = []
saouhsc= []
start = []
end = []
window = []
for line in open(sys.argv[1]):
    split=line.split('\t')
    if len(split) > 4:
        if split[3][0].isdigit() == True:
            inputreads.append(float(split[3]))
            saouhsc.append(split[0])
            start.append(split[1])
            end.append(split[2])
            window.append(split[4].rstrip('\n'))
            TActrl += 1
            Readsctrl += float(split[3])
            if float(split[3]) > 0: HitSitesctrl += 1

TAexp = 0; Readsexp = 0; HitSitesexp = 0
for line in open(sys.argv[2]):
    split=line.split()
    if len(split) > 3:
        if split[3][0].isdigit() == True:
            TAexp += 1
            Readsexp += float(split[3])
            if float(split[3]) > 0: HitSitesexp += 1

TAproportion =  float(HitSitesexp)/HitSitesctrl   

inputproportion = []
for i in inputreads:
	inputproportion.append(float(i)/Readsctrl)
inputproportiontanorm = []
for i in inputproportion:
	inputproportiontanorm.append(float(i)*TAproportion)

inputproportiontanorm.append(TAproportion-1)

multinominputsample = numpy.random.multinomial(Readsexp,inputproportiontanorm,100)
multinominputsample = multinominputsample.T
for i in range(len(multinominputsample[0])):
    summulti = sum(multinominputsample[:,i])
multisum = numpy.repeat(summulti,100)
difference=numpy.repeat(multisum-multinominputsample[-1,0],100)
correctionfactor = round(float(Readsexp)/difference[0],4)
correctedinput = [i*correctionfactor for i in multinominputsample]
bootstrapcontrol = numpy.delete(correctedinput, (-1), axis=0)
avgbootstrapcontrol = bootstrapcontrol.mean(axis=1)

x = 0
for x in range(0,len(avgbootstrapcontrol)):
    print saouhsc[x], '\t', start[x], '\t', end[x], '\t', avgbootstrapcontrol[x], '\t', window[x]
    x+=1


